#include <iostream>
#include "need-a-name.hh"
#include "subproject-base.hh"
#include "../algorithm/includes/algorithm.hh"
#include <string>
// #define DEBUG_GOTHROUGH

namespace GoThroughEntity{ 
  template<class S, class G, class M, class E, class I, class T>
  void interface(std::string flag, const S& subproject, const G& gridview, const M& mapper, const E& element, const I& i, const T& task, double elevolume, double* globalcoord, double rdt){
    const int dimworld = G::dimensionworld;
    int interscounter[gridsetting::dimension];
    int currentdim=0;    
    for(int tmpi=0;tmpi<gridsetting::dimension;tmpi++)
      interscounter[tmpi]=0;
    for(const auto& intersection : intersections(gridview,element)){
      auto intersectionGeometry = intersection.geometry();
      double intersectionArea = intersectionGeometry.volume();
      double pillarlen=elevolume/intersectionArea;
      Dune::FieldVector<double,dimworld> intersectionCenter = intersectionGeometry.center();
      const auto& intersectionReferenceElement = Dune::ReferenceElements<double,dimworld-1>::general(intersection.type());
      Dune::FieldVector<double,dimworld-1> intersectionLocalCenter = intersectionReferenceElement.position(0,0);
      Dune::FieldVector<double,dimworld> unitOuterNormal = intersection.unitOuterNormal(intersectionLocalCenter);
      double* unitouternormal = unitOuterNormal.data();
      typename M::Index j=0;
      for(int tmpi=0; tmpi<gridsetting::dimension; tmpi++)
	if(unitouternormal[tmpi]!=0)
	  currentdim=tmpi;
      interscounter[currentdim]+=1;

      //          std::cout << "enter interface: " << flag << std::endl;
#ifdef DEBUG_GOTHROUGH
      std::cout<< "interface currentdim: "<< currentdim<<std::endl;
      std::cout<< "interface currentdim direction: "<< unitouternormal[currentdim]<<std::endl;      
      std::cout<< "set initial condition-interface"<<std::endl;
#endif      
      if(memory::dgflag==1 && flag=="initialization"){
	// std::cout<< "dimworld: " << dimworld << std::endl;
	subproject->init->initial_function(globalcoord,pillarlen,currentdim,i);
	for(int tmpi=0; tmpi<std::max(1,memory::dgorder)*memory::phynumber+dataanalysis::auxnumber;tmpi++)
	  memory::variable_tmp[tmpi][i]=memory::variable[tmpi][i];
      }
      if(flag=="dt")
	task->get_dt(subproject,pillarlen,i);
      if(flag=="slope" || flag=="fvevolve" || flag=="dgevolve_calculateFlux"){
	if(intersection.neighbor())
	  j = mapper.index(intersection.outside());
	if(intersection.boundary())
	  j = -1;
      }
      if(flag=="slope")
	task->get_slope(subproject,globalcoord,unitouternormal,pillarlen,i,j,currentdim,interscounter);
      if(flag=="fvevolve"){
	// std::cout << "enter interface: " << flag << std::endl;      	
	task->fv_evolve(subproject,globalcoord,unitouternormal,pillarlen,i,j,currentdim,interscounter,rdt);
      }
      if(flag=="dgevolve_finiteElement")
	task->dg_evolve_finite_element(subproject,globalcoord,unitouternormal,pillarlen,i,j,currentdim,interscounter,rdt);
      if(flag=="dgevolve_calculateFlux")
	task->dg_evolve_calculate_flux(subproject,globalcoord,unitouternormal,pillarlen,i,j,currentdim,interscounter,rdt);
      if(flag=="sfeevolve")
        task->sfe_evolve(subproject,globalcoord,unitouternormal,pillarlen,i,j,currentdim,interscounter,rdt);  
    }
  }

  template<class S, class G, class M, class T>
  void element(std::string flag, const S& subproject, const G& gridview, const M& mapper, const T& task, double rdt){
    const int dimworld = G::dimensionworld;
    typedef typename G::ctype ct;
#ifdef DEBUG
#endif    
    std::cout << "enter element: " << flag << std::endl;
    for(const auto& element : Dune::elements(gridview,Dune::Partitions::all)){      
      const Dune::Geometry geo = element.geometry();
      double elevolume = geo.volume();
      typename M::Index i = mapper.index(element);
      Dune::FieldVector<ct,dimworld> global = geo.center();
      double globalcoord[dimworld];
      for(int tmpi=0; tmpi<dimworld; tmpi++)
	globalcoord[tmpi] = global[tmpi];
#ifdef DEBUG_GOTHROUGH
      std::cout << "element index: " << i << std::endl;	
#endif	
      // tasks only need element-loop
      if(memory::dgflag==0 && flag=="initialization"){
	double dummyPillarlen=1.0;
	int dummyCurrentlen=0;	
	subproject->init->initial_function(globalcoord,dummyPillarlen,dummyCurrentlen,i);
	for(int tmpi=0; tmpi<std::max(1,memory::dgorder)*memory::phynumber+dataanalysis::auxnumber;tmpi++)
	  memory::variable_tmp[tmpi][i]=memory::variable[tmpi][i];
	task->dummytask();
      }
      if(memory::dgflag==1 && flag=="dgevolve_combineWeight"){
		task->dg_evolve_combine_weight(i);
      }
      // tasks also need intersection-loop
      if(memory::dgflag==1 && flag=="initialization"){
        // std::cout<< "set initial condition-element"<<std::endl;	
	interface(flag,subproject,gridview,mapper,element,i,task,elevolume,globalcoord,rdt);	
      }
      if(flag=="dt" || flag=="slope" || flag=="fvevolve" || flag=="dgevolve_finiteElement" || flag=="dgevolve_calculateFlux"|| flag=="sfeevolve"){
	interface(flag,subproject,gridview,mapper,element,i,task,elevolume,globalcoord,rdt);
      }
    }
  }
  
  ///////////////////////////////////////////////////////
  /////////////////   TIME ADVANCE    ///////////////////
  ///////////////////////////////////////////////////////
  
  template<class S, class G, class M>
  void exchange_datum(std::string flag, const S& subproject, const G& gridview, const M& mapper, std::vector<std::vector<double> >& exchange_var){
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
      VectorExchange<M,std::vector<double> > v(mapper,exchange_var[tmpi]);
      gridview.template communicate<VectorExchange<M,std::vector<double> > >(v,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
    }
  }
  template<class S, class G, class M>
  void exchange_datum2(std::string flag, const S& subproject, const G& gridview, const M& mapper, std::vector<std::vector<std::vector<double> > >& exchange_var){
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
      for(int tmpj=0; tmpj<gridsetting::dimension; tmpj++){
	VectorExchange<M,std::vector<double> > v(mapper,exchange_var[tmpi][tmpj]);
	gridview.template communicate<VectorExchange<M,std::vector<double> > >(v,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
      }
    }
  }

  template<class S, class G, class M, class T>
  void timeadvance_euler(std::string flag, const S& subproject, const G& gridview, const M& mapper, const T& task, double dt){
    std::cout << "euler time advance" << std::endl;
    // std::cout << "task flag: " << flag  << std::endl;    
    timerule::roundnumber=0;
    element(flag,subproject,gridview,mapper,task,dt);
    for(int tmpj=0; tmpj<mapper.size(); ++tmpj)
      for(int tmpd=0; tmpd<std::max(1,memory::dgorder); tmpd++)
	for(int tmpi=0; tmpi<memory::phynumber; tmpi++){	
	  memory::variable_tmp[tmpd*memory::phynumber+tmpi][tmpj]=memory::variable[tmpd*memory::phynumber+tmpi][tmpj]+dt*memory::update[timerule::roundnumber][tmpd*memory::phynumber+tmpi][tmpj];
	  memory::variable[tmpd*memory::phynumber+tmpi][tmpj]=memory::variable_tmp[tmpd*memory::phynumber+tmpi][tmpj];
      }
    // exchange_datum(flag,subproject,gridview,mapper,memory::variable);
    exchange_datum(flag,subproject,gridview,mapper,memory::variable_tmp);
  };

  template<class S, class G, class M, class T>
  void timeadvance_rk2(std::string flag, const S& subproject, const G& gridview, const M& mapper, const T& task, double dt){
    std::cout << "2-order runge-kutta -- predictor-corrector" << std::endl;
    // k1
    timerule::roundnumber=0;
    element(flag,subproject,gridview,mapper,task,dt);
    // k2
    for(int tmpj=0; tmpj<mapper.size(); ++tmpj)
      for(int tmpd=0; tmpd<std::max(1,memory::dgorder); tmpd++)
	for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
	  memory::variable_tmp[tmpd*memory::phynumber+tmpi][tmpj]=memory::variable[tmpd*memory::phynumber+tmpi][tmpj]+dt*memory::update[timerule::roundnumber][tmpd*memory::phynumber+tmpi][tmpj];
    exchange_datum(flag,subproject,gridview,mapper,memory::variable_tmp);    
    timerule::roundnumber=1;
    element(flag,subproject,gridview,mapper,task,dt);
    // dt/2*(k1+k2)
    for(int tmpj=0; tmpj<mapper.size(); ++tmpj)
      for(int tmpd=0; tmpd<std::max(1,memory::dgorder); tmpd++)
	for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
	  memory::variable_tmp[tmpd*memory::phynumber+tmpi][tmpj]=memory::variable[tmpd*memory::phynumber+tmpi][tmpj]+dt/2.0*(memory::update[timerule::roundnumber-1][tmpd*memory::phynumber+tmpi][tmpj]+memory::update[timerule::roundnumber][tmpd*memory::phynumber+tmpi][tmpj]);
	  memory::variable[tmpd*memory::phynumber+tmpi][tmpj]=memory::variable_tmp[tmpd*memory::phynumber+tmpi][tmpj];
      }
    // exchange_datum(flag,subproject,gridview,mapper,memory::variable);
    exchange_datum(flag,subproject,gridview,mapper,memory::variable_tmp);    
  };

  template<class S, class G, class M, class T>
  void timeadvance_rk3(std::string flag, const S& subproject, const G& gridview, const M& mapper, const T& task, double dt){
    std::cout << "3-order runge-kutta -- Heun" << std::endl;
    // k1
    timerule::roundnumber=0;
    element(flag,subproject,gridview,mapper,task,dt);
    // k2
    for(int tmpj=0; tmpj<mapper.size(); ++tmpj)
      for(int tmpd=0; tmpd<std::max(1,memory::dgorder); tmpd++)
	for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
	  memory::variable_tmp[tmpd*memory::phynumber+tmpi][tmpj]=memory::variable[tmpd*memory::phynumber+tmpi][tmpj]+dt/3.0*memory::update[timerule::roundnumber][tmpd*memory::phynumber+tmpi][tmpj];
    exchange_datum(flag,subproject,gridview,mapper,memory::variable_tmp); 
    timerule::roundnumber=1;
    element(flag,subproject,gridview,mapper,task,dt);
    // k3
    for(int tmpj=0; tmpj<mapper.size(); ++tmpj)
      for(int tmpd=0; tmpd<std::max(1,memory::dgorder); tmpd++)
	for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
	  memory::variable_tmp[tmpd*memory::phynumber+tmpi][tmpj]=memory::variable[tmpd*memory::phynumber+tmpi][tmpj]+dt*2.0/3.0*memory::update[timerule::roundnumber][tmpd*memory::phynumber+tmpi][tmpj];
    exchange_datum(flag,subproject,gridview,mapper,memory::variable_tmp); 
    timerule::roundnumber=2;
    element(flag,subproject,gridview,mapper,task,dt);
    // dt/4*(k1+3*k3)
    for(int tmpj=0; tmpj<mapper.size(); ++tmpj)
      for(int tmpd=0; tmpd<std::max(1,memory::dgorder); tmpd++)
	for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
	  memory::variable_tmp[tmpd*memory::phynumber+tmpi][tmpj]=memory::variable[tmpd*memory::phynumber+tmpi][tmpj]+dt/4.0*(memory::update[0][tmpd*memory::phynumber+tmpi][tmpj]+3.0*memory::update[2][tmpd*memory::phynumber+tmpi][tmpj]);
	  memory::variable[tmpd*memory::phynumber+tmpi][tmpj]=memory::variable_tmp[tmpd*memory::phynumber+tmpi][tmpj];
      }
    exchange_datum(flag,subproject,gridview,mapper,memory::variable_tmp);    
  };

  template<class S, class G, class M, class T>
  void timeadvance_rk4(std::string flag, const S& subproject, const G& gridview, const M& mapper, const T& task, double dt){
    std::cout << "4-order classic runge-kutta" << std::endl;
    element(flag,subproject,gridview,mapper,task,dt);    
  };

  //  template<class S, class G, class M>  
  // void weight_to(std::string flag, const S& subproject, const G& gridview, const M& mapper){
    // element(flag,subproject,gridview,mapper,task,dt);        
    //    double basis_middle[memory::dgorder];
    // basis.getBasisMiddle(basis_middle);
    //for(int tmpj=0; tmpj<mapper.size(); ++tmpj){
    // }
    //   for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
    // 	// memory::variable_tmp[(memory::dgorder-1)*memory::phynumber+tmpi][tmpj]=0.0;
    // 	//memory::variable[tmpi][tmpj]=0.0;	
    // 	//for(int tmpd=0; tmpd<memory::dgorder; tmpd++){
    // 	//  memory::variable_tmp[(memory::dgorder-1)*memory::phynumber+tmpi][tmpj]+=memory::variable_tmp[tmpd*memory::phynumber+tmpi][tmpj]*basis_middle[tmpd];
    // 	//}
    //     //memory::variable[(memory::dgorder-1)*memory::phynumber+tmpi][tmpj]=memory::variable_tmp[(memory::dgorder-1)*memory::phynumber+tmpi][tmpj];	
    //   }
  //};
  ///////////////////////////////////////////////////////
  ////////////////////    GO    /////////////////////////
  ///////////////////////////////////////////////////////
  
  template<class S, class G, class M>
  void go(std::string flag, const S& subproject, const G& gridview, const M& mapper){
    Dummy* dummy=new Dummy();
    TimeStep* timestep=new TimeStep();
    Slope* slope=new Slope();    
    FVEvolve* fvevolve=new FVEvolve();
    DGEvolveFiniteElement* dgevolvefiniteelement=new DGEvolveFiniteElement();    
    DGEvolveCalculateFlux* dgevolvecalculateflux=new DGEvolveCalculateFlux();
    DGEvolveCombineWeight* dgevolvecombineweight=new DGEvolveCombineWeight();    
    SFEEvolve* sfeevolve=new SFEEvolve();
    if(flag=="initialization"){      
      std::cout<< "set initial condition"  <<std::endl;
      element(flag,subproject,gridview,mapper,dummy,0.0);
    }
    if(flag=="dt") {
      std::cout<< "calculate time step dt"  <<std::endl;
      element(flag,subproject,gridview,mapper,timestep,0.0);
      //TODO: dt needs exchange across processes      
      timerule::dt = timestep->global_dt;
#ifdef DEBUG
      std::cout << "timerule::dt: " << timerule::dt << std::endl;
#endif
    }
    if(flag=="fvevolve" || flag=="dgevolve_finiteElement" || flag=="dgevolve_calculateFlux" || flag=="sfeevolve"){
      for(int tmpi=0; tmpi<memory::updatenumber; tmpi++)
	for(int tmpj=0; tmpj<std::max(1,memory::dgorder)*memory::phynumber; tmpj++)
	  std::fill(memory::update[tmpi][tmpj].begin(),memory::update[tmpi][tmpj].end(),0.0);      
    }
    if(flag=="slope") {
      element(flag,subproject,gridview,mapper,slope,0.0);
      exchange_datum2(flag,subproject,gridview,mapper,memory::limit_slope);
    }
    if(flag=="fvevolve") {
      std::cout<< "fv time advance"  <<std::endl;
      switch(setting::timeadvance){
      case setting::euler:
	timeadvance_euler(flag,subproject,gridview,mapper,fvevolve,timerule::dt);
	break;
      case setting::rk2:
	timeadvance_rk2(flag,subproject,gridview,mapper,fvevolve,timerule::dt);
	break;
      case setting::rk3:
	timeadvance_rk3(flag,subproject,gridview,mapper,fvevolve,timerule::dt);
	break;
      }
    }

    if(flag=="dgevolve_finiteElement") {
      std::cout<< "dg time advance for finite element"  <<std::endl;
      switch(setting::timeadvance){
      case setting::euler:
	timeadvance_euler(flag,subproject,gridview,mapper,dgevolvefiniteelement,timerule::dt);
	break;
      case setting::rk2:
	timeadvance_rk2(flag,subproject,gridview,mapper,dgevolvefiniteelement,timerule::dt);
	break;
      case setting::rk3:
	timeadvance_rk3(flag,subproject,gridview,mapper,dgevolvefiniteelement,timerule::dt);
	break;
      }
    }    
    
    if(flag=="dgevolve_calculateFlux") {
      std::cout<< "dg time advance for flux calculation"  <<std::endl;
      switch(setting::timeadvance){
      case setting::euler:
	timeadvance_euler(flag,subproject,gridview,mapper,dgevolvecalculateflux,timerule::dt);
	break;
      case setting::rk2:
	timeadvance_rk2(flag,subproject,gridview,mapper,dgevolvecalculateflux,timerule::dt);
	break;
      case setting::rk3:
	timeadvance_rk3(flag,subproject,gridview,mapper,dgevolvecalculateflux,timerule::dt);
	break;
      }
    }

    if(flag=="dgevolve_combineWeight") {
      element(flag,subproject,gridview,mapper,dgevolvecombineweight,0.0);      
    }

    
    // if(flag=="sfeevolve") {
    // }
  }
  
};
